Genome-scale metabolic modeling of the human milk oligosaccharide utilization by Bifidobacterium longum subsp. infantis.

Bifidobacterium longum subsp. infantis is a representative and dominant species in the infant gut and is considered a beneficial microbe. This organism displays multiple adaptations to thrive in the infant gut, regarded as a model for human milk oligosaccharides (HMOs) utilization. These carbohydrates are abundant in breast milk and include different molecules based on lactose. They contain fucose, sialic acid, and N-acetylglucosamine. Bifidobacterium metabolism is complex, and a systems view of relevant metabolic pathways and exchange metabolites during HMO consumption is missing. To address this limitation, a refined genome-scale network reconstruction of this bacterium is presented using a previous reconstruction of B. infantis ATCC 15967 as a template. The latter was expanded based on an extensive revision of genome annotations, current literature, and transcriptomic data integration. The metabolic reconstruction (iLR578) accounted for 578 genes, 1,047 reactions, and 924 metabolites. Starting from this reconstruction, we built context-specific genome-scale metabolic models using RNA-seq data from cultures growing in lactose and three HMOs. The models revealed notable differences in HMO metabolism depending on the functional characteristics of the substrates. Particularly, fucosyl-lactose showed a divergent metabolism due to a fucose moiety. High yields of lactate and acetate were predicted under growth rate maximization in all conditions, whereas formate, ethanol, and 1,2-propanediol were substantially lower. Similar results were also obtained under near-optimal growth on each substrate when varying the empirically observed acetate-to-lactate production ratio. Model predictions displayed reasonable agreement between central carbon metabolism fluxes and expression data across all conditions. Flux coupling analysis revealed additional connections between succinate exchange and arginine and sulfate metabolism and a strong coupling between central carbon reactions and adenine metabolism. More importantly, specific networks of coupled reactions under each carbon source were derived and analyzed. Overall, the presented network reconstruction constitutes a valuable platform for probing the metabolism of this prominent infant gut bifidobacteria.IMPORTANCEThis work presents a detailed reconstruction of the metabolism of Bifidobacterium longum subsp. infantis, a prominent member of the infant gut microbiome, providing a systems view of its metabolism of human milk oligosaccharides.

LooplessFluxSampler: an efficient toolbox for sampling the loopless flux solution space of metabolic models.

BACKGROUND: Uniform random sampling of mass-balanced flux solutions offers an unbiased appraisal of the capabilities of metabolic networks. Unfortunately, it is impossible to avoid thermodynamically infeasible loops in flux samples when using convex samplers on large metabolic models. Current strategies for randomly sampling the non-convex loopless flux space display limited efficiency and lack theoretical guarantees. RESULTS: Here, we present LooplessFluxSampler, an efficient algorithm for exploring the loopless mass-balanced flux solution space of metabolic models, based on an Adaptive Directions Sampling on a Box (ADSB) algorithm. ADSB is rooted in the general Adaptive Direction Sampling (ADS) framework, specifically the Parallel ADS, for which theoretical convergence and irreducibility results are available for sampling from arbitrary distributions. By sampling directions that adapt to the target distribution, ADSB traverses more efficiently the sample space achieving faster mixing than other methods. Importantly, the presented algorithm is guaranteed to target the uniform distribution over convex regions, and it provably converges on the latter distribution over more general (non-convex) regions provided the sample can have full support. CONCLUSIONS: LooplessFluxSampler enables scalable statistical inference of the loopless mass-balanced solution space of large metabolic models. Grounded in a theoretically sound framework, this toolbox provides not only efficient but also reliable results for exploring the properties of the almost surely non-convex loopless flux space. Finally, LooplessFluxSampler includes a Markov Chain diagnostics suite for assessing the quality of the final sample and the performance of the algorithm.

Microbial interactions and the homeostasis of the gut microbiome: the role of Bifidobacterium.

The human gut is home to trillions of microorganisms that influence several aspects of our health. This dense microbial community targets almost all dietary polysaccharides and releases multiple metabolites, some of which have physiological effects on the host. A healthy equilibrium between members of the gut microbiota, its microbial diversity, and their metabolites is required for intestinal health, promoting regulatory or anti-inflammatory immune responses. In contrast, the loss of this equilibrium due to antibiotics, low fiber intake, or other conditions results in alterations in gut microbiota composition, a term known as gut dysbiosis. This dysbiosis can be characterized by a reduction in health-associated microorganisms, such as butyrate-producing bacteria, enrichment of a small number of opportunistic pathogens, or a reduction in microbial diversity. Bifidobacterium species are key species in the gut microbiome, serving as primary degraders and contributing to a balanced gut environment in various ways. Colonization resistance is a fundamental property of gut microbiota for the prevention and control of infections. This community competes strongly with foreign microorganisms, such as gastrointestinal pathogens, antibiotic-resistant bacteria, or even probiotics. Resistance to colonization is based on microbial interactions such as metabolic cross-feeding, competition for nutrients, or antimicrobial-based inhibition. These interactions are mediated by metabolites and metabolic pathways, representing the inner workings of the gut microbiota, and play a protective role through colonization resistance. This review presents a rationale for how microbial interactions provide resistance to colonization and gut dysbiosis, highlighting the protective role of Bifidobacterium species.

The topology of genome-scale metabolic reconstructions unravels independent modules and high network flexibility.

The topology of metabolic networks is recognisably modular with modules weakly connected apart from sharing a pool of currency metabolites. Here, we defined modules as sets of reversible reactions isolated from the rest of metabolism by irreversible reactions except for the exchange of currency metabolites. Our approach identifies topologically independent modules under specific conditions associated with different metabolic functions. As case studies, the E.coli iJO1366 and Human Recon 2.2 genome-scale metabolic models were split in 103 and 321 modules respectively, displaying significant correlation patterns in expression data. Finally, we addressed a fundamental question about the metabolic flexibility conferred by reversible reactions: "Of all Directed Topologies (DTs) defined by fixing directions to all reversible reactions, how many are capable of carrying flux through all reactions?". Enumeration of the DTs for iJO1366 model was performed using an efficient depth-first search algorithm, rejecting infeasible DTs based on mass-imbalanced and loopy flux patterns. We found the direction of 79% of reversible reactions must be defined before all directions in the network can be fixed, granting a high degree of flexibility.

Metabolic Modeling of Wine Fermentation at Genome Scale.

Wine fermentation is an ancient biotechnological process mediated by different microorganisms such as yeast and bacteria. Understanding of the metabolic and physiological phenomena taking place during this process can be now attained at a genome scale with the help of metabolic models. In this chapter, we present a detailed protocol for modeling wine fermentation using genome-scale metabolic models. In particular, we illustrate how metabolic fluxes can be computed, optimized and interpreted, for both yeast and bacteria under winemaking conditions. We also show how nutritional requirements can be determined and simulated using these models in relevant test cases. This chapter introduces fundamental concepts and practical steps for applying flux balance analysis in wine fermentation, and as such, it is intended for a broad microbiology audience as well as for practitioners in the metabolic modeling field.

GRASP: a computational platform for building kinetic models of cellular metabolism.

Summary: Kinetic models of metabolism are crucial to understand the inner workings of cell metabolism. By taking into account enzyme regulation, detailed kinetic models can provide accurate predictions of metabolic fluxes. Comprehensive consideration of kinetic regulation requires highly parameterized non-linear models, which are challenging to build and fit using available modelling tools. Here, we present a computational package implementing the GRASP framework for building detailed kinetic models of cellular metabolism. By defining the mechanisms of enzyme regulation and a reference state described by reaction fluxes and their corresponding Gibbs free energy ranges, GRASP can efficiently sample an arbitrarily large population of thermodynamically feasible kinetic models. If additional experimental data are available (fluxes, enzyme and metabolite concentrations), these can be integrated to generate models that closely reproduce these observations using an approximate Bayesian computation fitting framework. Within the same framework, model selection tasks can be readily performed. Availability and implementation: GRASP is implemented as an open-source package in the MATLAB environment. The software runs in Windows, macOS and Linux, is documented (graspk.readthedocs.io) and unit-tested. GRASP is freely available at github.com/biosustain/GRASP. Supplementary information: Supplementary data are available at Bioinformatics Advances online.

Screening of COVID-19 cases through a Bayesian network symptoms model and psychophysical olfactory test.

The sudden loss of smell is among the earliest and most prevalent symptoms of COVID-19 when measured with a clinical psychophysical test. Research has shown the potential impact of frequent screening for olfactory dysfunction, but existing tests are expensive and time consuming. We developed a low-cost ($0.50/test) rapid psychophysical olfactory test (KOR) for frequent testing and a model-based COVID-19 screening framework using a Bayes Network symptoms model. We trained and validated the model on two samples: suspected COVID-19 cases in five healthcare centers (n = 926; 33% prevalence, 309 RT-PCR confirmed) and healthy miners (n = 1,365; 1.1% prevalence, 15 RT-PCR confirmed). The model predicted COVID-19 status with 76% and 96% accuracy in the healthcare and miners samples, respectively (healthcare: AUC = 0.79 [0.75-0.82], sensitivity: 59%, specificity: 87%; miners: AUC = 0.71 [0.63-0.79], sensitivity: 40%, specificity: 97%, at 0.50 infection probability threshold). Our results highlight the potential for low-cost, frequent, accessible, routine COVID-19 testing to support society's reopening.

Inferring composition and function of the human gut microbiome in time and space: A review of genome-scale metabolic modelling tools.

The human gut hosts a complex community of microorganisms that directly influences gastrointestinal physiology, playing a central role in human health. Because of its importance, the metabolic interplay between the gut microbiome and host metabolism has gained special interest. While there has been great progress in the field driven by metagenomics and experimental studies, the mechanisms underpinning microbial composition and interactions in the microbiome remain poorly understood. Genome-scale metabolic models are mathematical structures capable of describing the metabolic potential of microbial cells. They are thus suitable tools for probing the metabolic properties of microbial communities. In this review, we discuss the most recent and relevant genome-scale metabolic modelling tools for inferring the composition, interactions, and ultimately, biological function of the constituent species of a microbial community with special emphasis in the gut microbiota. Particular attention is given to constraint-based metabolic modelling methods as well as hybrid agent-based methods for capturing the interactions and behavior of the community in time and space. Finally, we discuss the challenges hindering comprehensive modelling of complex microbial communities and its application for the in-silico design of microbial consortia with therapeutic functions.

Engineering Saccharomyces cerevisiae for the Overproduction of ß-Ionone and Its Precursor ß-Carotene.

ß-ionone is a commercially attractive industrial fragrance produced naturally from the cleavage of the pigment ß-carotene in plants. While the production of this ionone is typically performed using chemical synthesis, environmentally friendly and consumer-oriented biotechnological production is gaining increasing attention. A convenient cell factory to address this demand is the yeast Saccharomyces cerevisiae. However, current ß-ionone titers and yields are insufficient for commercial bioproduction. In this work, we optimized S. cerevisiae for the accumulation of high amounts of ß-carotene and its subsequent conversion to ß-ionone. For this task, we integrated systematically the heterologous carotenogenic genes (CrtE, CrtYB and CrtI) from Xanthophyllomyces dendrorhous using markerless genome editing CRISPR/Cas9 technology; and evaluated the transcriptional unit architecture (bidirectional or tandem), integration site, and impact of gene dosage, first on ß-carotene accumulation, and later, on ß-ionone production. A single-copy insertion of the carotenogenic genes in high expression loci of the wild-type yeast CEN.Pk2 strain yielded 4 mg/gDCW of total carotenoids, regardless of the transcriptional unit architecture employed. Subsequent fine-tuning of the carotenogenic gene expression enabled reaching 16 mg/gDCW of total carotenoids, which was further increased to 32 mg/gDCW by alleviating the known pathway bottleneck catalyzed by the hydroxymethylglutaryl-CoA reductase (HMGR1). The latter yield represents the highest total carotenoid concentration reported to date in S. cerevisiae for a constitutive expression system. For ß-ionone synthesis, single and multiple copies of the carotene cleavage dioxygenase 1 (CCD1) gene from Petunia hybrida (PhCCD1) fused with a membrane destination peptide were expressed in the highest ß-carotene-producing strains, reaching up to 33 mg/L of ß-ionone in the culture medium after 72-h cultivation in shake flasks. Finally, interrogation of a contextualized genome-scale metabolic model of the producer strains pointed to PhCCD1 unspecific cleavage activity as a potentially limiting factor reducing ß-ionone production. Overall, the results of this work constitute a step toward the industrial production of this ionone and, more broadly, they demonstrate that biotechnological production of apocarotenoids is technically feasible.

Genomic integration of unclonable gene expression cassettes in Saccharomyces cerevisiae using rapid cloning-free workflows.

Most DNA assembly methods require bacterial amplification steps, which restrict its application to genes that can be cloned in the bacterial host without significant toxic effects. However, genes that cannot be cloned in bacteria do not necessarily exert toxic effects on the final host. In order to tackle this issue, we adapted two DNA assembly workflows for rapid, cloning-free construction and genomic integration of expression cassettes in Saccharomyces cerevisiae. One method is based on a modified Gibson assembly, while the other relies on a direct assembly and integration of linear PCR products by yeast homologous recombination. The methods require few simple experimental steps, and their performance was evaluated for the assembly and integration of unclonable zeaxanthin epoxidase expression cassettes in yeast. Results showed that up to 95% integration efficiency can be reached with minimal experimental effort. The presented workflows can be employed as rapid gene integration tools for yeast, especially tailored for integrating unclonable genes.

Contextualized genome-scale model unveils high-order metabolic effects of the specific growth rate and oxygenation level in recombinant Pichia pastoris.

Pichia pastoris is recognized as a biotechnological workhorse for recombinant protein expression. The metabolic performance of this microorganism depends on genetic makeup and culture conditions, amongst which the specific growth rate and oxygenation level are critical. Despite their importance, only their individual effects have been assessed so far, and thus their combined effects and metabolic consequences still remain to be elucidated. In this work, we present a comprehensive framework for revealing high-order (i.e., individual and combined) metabolic effects of the above parameters in glucose-limited continuous cultures of P. pastoris, using thaumatin production as a case study. Specifically, we employed a rational experimental design to calculate statistically significant metabolic effects from multiple chemostat data, which were later contextualized using a refined and highly predictive genome-scale metabolic model of this yeast under the simulated conditions. Our results revealed a negative effect of the oxygenation on the specific product formation rate (thaumatin), and a positive effect on the biomass yield. Notably, we identified a novel positive combined effect of both the specific growth rate and oxygenation level on the specific product formation rate. Finally, model predictions indicated an opposite relationship between the oxygenation level and the growth-associated maintenance energy (GAME) requirement, suggesting a linear GAME decrease of 0.56 mmol ATP/gDCW per each 1% increase in oxygenation level, which translated into a 44% higher metabolic cost under low oxygenation compared to high oxygenation. Overall, this work provides a systematic framework for mapping high-order metabolic effects of different culture parameters on the performance of a microbial cell factory. Particularly in this case, it provided valuable insights about optimal operational conditions for protein production in P. pastoris.

Build Your Bioprocess on a Solid Strain-ß-Carotene Production in Recombinant Saccharomyces cerevisiae.

Robust fermentation performance of microbial cell factories is critical for successful scaling of a biotechnological process. From shake flask cultivations to industrial-scale bioreactors, consistent strain behavior is fundamental to achieve the production targets. To assert the importance of this feature, we evaluated the impact of the yeast strain design and construction method on process scalability -from shake flasks to bench-scale fed-batch fermentations- using two recombinant Saccharomyces cerevisiae strains capable of producing ß-carotene; SM14 and ßcar1.2 strains. SM14 strain, obtained previously from adaptive evolution experiments, was capable to accumulate up to 21 mg/gDCW of ß-carotene in 72 h shake flask cultures; while the ßcar1.2, constructed by overexpression of carotenogenic genes, only accumulated 5.8 mg/gDCW of carotene. Surprisingly, fed-batch cultivation of these strains in 1L bioreactors resulted in opposite performances. ßcar1.2 strain reached much higher biomass and ß-carotene productivities (1.57 g/L/h and 10.9 mg/L/h, respectively) than SM14 strain (0.48 g/L/h and 3.1 mg/L/h, respectively). Final ß-carotene titers were 210 and 750 mg/L after 80 h cultivation for SM14 and ßcar1.2 strains, respectively. Our results indicate that these substantial differences in fermentation parameters are mainly a consequence of the exacerbated Crabtree effect of the SM14 strain. We also found that the strategy used to integrate the carotenogenic genes into the chromosomes affected the genetic stability of strains, although the impact was significantly minor. Overall, our results indicate that shake flasks fermentation parameters are poor predictors of the fermentation performance under industrial-like conditions, and that appropriate construction designs and performance tests must be conducted to properly assess the scalability of the strain and the bioprocess.

Expanding Metabolic Capabilities Using Novel Pathway Designs: Computational Tools and Case Studies.

Design and selection of efficient metabolic pathways is critical for the success of metabolic engineering endeavors. Convenient pathways should not only produce the target metabolite in high yields but also are required to be thermodynamically feasible under production conditions, and to prefer efficient enzymes. To support the design and selection of such pathways, different computational approaches have been proposed for exploring the feasible pathway space under many of the above constraints. In this review, an overview of recent constraint-based optimization frameworks for metabolic pathway prediction, as well as relevant pathway engineering case studies that highlight the importance of rational metabolic designs is presented. Despite the availability and suitability of in silico design tools for metabolic pathway engineering, scarce-although increasing-application of computational outcomes is found. Finally, challenges and limitations hindering the broad adoption and successful application of these tools in metabolic engineering projects are discussed.

From reconstruction to C4 metabolic engineering: A case study for overproduction of polyhydroxybutyrate in bioenergy grasses.

The compartmentalization of C4 plants increases photosynthetic efficiency, while constraining how material and energy must flow in leaf tissues. To capture this metabolic phenomenon, a generic plant metabolic reconstruction was replicated into four connected spatiotemporal compartments, namely bundle sheath (B) and mesophyll (M) across the day and night cycle. The C4 leaf model was used to explore how amenable polyhydroxybutyrate (PHB) production is with these four compartments working cooperatively. A strategic pattern of metabolite conversion and exchange emerged from a systems-level network that has very few constraints imposed; mainly the sequential two-step carbon capture in mesophyll, then bundle sheath and photosynthesis during the day only. The building of starch reserves during the day and their mobilization during the night connects day and night metabolism. Flux simulations revealed that PHB production did not require rerouting of metabolic pathways beyond what is already utilised for growth. PHB yield was sensitive to photoassimilation capacity, availability of carbon reserves, ATP maintenance, relative photosynthetic activity of B and M, and type of metabolites exchanged in the plasmodesmata, but not sensitive towards compartmentalization. Hence, the compartmentalization issues currently encountered are likely to be kinetic or thermodynamic limitations rather than stoichiometric.

Formulation, construction and analysis of kinetic models of metabolism: A review of modelling frameworks.

Kinetic models are critical to predict the dynamic behaviour of metabolic networks. Mechanistic kinetic models for large networks remain uncommon due to the difficulty of fitting their parameters. Recent modelling frameworks promise new ways to overcome this obstacle while retaining predictive capabilities. In this review, we present an overview of the relevant mathematical frameworks for kinetic formulation, construction and analysis. Starting with kinetic formalisms, we next review statistical methods for parameter inference, as well as recent computational frameworks applied to the construction and analysis of kinetic models. Finally, we discuss opportunities and limitations hindering the development of larger kinetic reconstructions.

Fast-SNP: a fast matrix pre-processing algorithm for efficient loopless flux optimization of metabolic models.

MOTIVATION: Computation of steady-state flux solutions in large metabolic models is routinely performed using flux balance analysis based on a simple LP (Linear Programming) formulation. A minimal requirement for thermodynamic feasibility of the flux solution is the absence of internal loops, which are enforced using 'loopless constraints'. The resulting loopless flux problem is a substantially harder MILP (Mixed Integer Linear Programming) problem, which is computationally expensive for large metabolic models. RESULTS: We developed a pre-processing algorithm that significantly reduces the size of the original loopless problem into an easier and equivalent MILP problem. The pre-processing step employs a fast matrix sparsification algorithm-Fast- sparse null-space pursuit (SNP)-inspired by recent results on SNP. By finding a reduced feasible 'loop-law' matrix subject to known directionalities, Fast-SNP considerably improves the computational efficiency in several metabolic models running different loopless optimization problems. Furthermore, analysis of the topology encoded in the reduced loop matrix enabled identification of key directional constraints for the potential permanent elimination of infeasible loops in the underlying model. Overall, Fast-SNP is an effective and simple algorithm for efficient formulation of loop-law constraints, making loopless flux optimization feasible and numerically tractable at large scale. AVAILABILITY AND IMPLEMENTATION: Source code for MATLAB including examples is freely available for download at http://www.aibn.uq.edu.au/cssb-resources under Software. Optimization uses Gurobi, CPLEX or GLPK (the latter is included with the algorithm). CONTACT: lars.nielsen@uq.edu.auSupplementary information: Supplementary data are available at Bioinformatics online.

Construction of feasible and accurate kinetic models of metabolism: A Bayesian approach.

Kinetic models are essential to quantitatively understand and predict the behaviour of metabolic networks. Detailed and thermodynamically feasible kinetic models of metabolism are inherently difficult to formulate and fit. They have a large number of heterogeneous parameters, are non-linear and have complex interactions. Many powerful fitting strategies are ruled out by the intractability of the likelihood function. Here, we have developed a computational framework capable of fitting feasible and accurate kinetic models using Approximate Bayesian Computation. This framework readily supports advanced modelling features such as model selection and model-based experimental design. We illustrate this approach on the tightly-regulated mammalian methionine cycle. Sampling from the posterior distribution, the proposed framework generated thermodynamically feasible parameter samples that converged on the true values, and displayed remarkable prediction accuracy in several validation tests. Furthermore, a posteriori analysis of the parameter distributions enabled appraisal of the systems properties of the network (e.g., control structure) and key metabolic regulations. Finally, the framework was used to predict missing allosteric interactions.

ll-ACHRB: a scalable algorithm for sampling the feasible solution space of metabolic networks.

MOTIVATION: Random sampling of the solution space has emerged as a popular tool to explore and infer properties of large metabolic networks. However, conventional sampling approaches commonly used do not eliminate thermodynamically unfeasible loops. RESULTS: In order to overcome this limitation, we developed an efficient sampling algorithm called loopless Artificially Centered Hit-and-Run on a Box (ll-ACHRB). This algorithm is inspired by the Hit-and-Run on a Box algorithm for uniform sampling from general regions, but employs the directions of choice approach of Artificially Centered Hit-and-Run. A novel strategy for generating feasible warmup points improved both sampling efficiency and mixing. ll-ACHRB shows overall better performance than current strategies to generate feasible flux samples across several models. Furthermore, we demonstrate that a failure to eliminate unfeasible loops greatly affects sample statistics, in particular the correlation structure. Finally, we discuss recommendations for the interpretation of sampling results and possible algorithmic improvements. AVAILABILITY AND IMPLEMENTATION: Source code for MATLAB and OCTAVE including examples are freely available for download at http://www.aibn.uq.edu.au/cssb-resources under Software. Optimization runs can use Gurobi Optimizer (by default if available) or GLPK (included with the algorithm). CONTACT: lars.nielsen@uq.edu.au SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

A probabilistic framework for the exploration of enzymatic capabilities based on feasible kinetics and control analysis.

BACKGROUND: Analysis of limiting steps within enzyme-catalyzed reactions is fundamental to understand their behavior and regulation. Methods capable of unravelling control properties and exploring kinetic capabilities of enzymatic reactions would be particularly useful for protein and metabolic engineering. While single-enzyme control analysis formalism has previously been applied to well-studied enzymatic mechanisms, broader application of this formalism is limited in practice by the limited amount of kinetic data and the difficulty of describing complex allosteric mechanisms. METHODS: To overcome these limitations, we present here a probabilistic framework enabling control analysis of previously unexplored mechanisms under uncertainty. By combining a thermodynamically consistent parameterization with an efficient Sequential Monte Carlo sampler embedded in a Bayesian setting, this framework yields insights into the capabilities of enzyme-catalyzed reactions with modest kinetic information, provided that the catalytic mechanism and a thermodynamic reference point are defined. RESULTS: The framework was used to unravel the impact of thermodynamic affinity, substrate saturation levels and effector concentrations on the flux control and response coefficients of a diverse set of enzymatic reactions. CONCLUSIONS: Our results highlight the importance of the metabolic context in the control analysis of isolated enzymes as well as the use of statistically sound methods for their interpretation. GENERAL SIGNIFICANCE: This framework significantly expands our current capabilities for unravelling the control properties of general reaction kinetics with limited amount of information. This framework will be useful for both theoreticians and experimentalists in the field.

A multi-tissue genome-scale metabolic modeling framework for the analysis of whole plant systems.

Genome scale metabolic modeling has traditionally been used to explore metabolism of individual cells or tissues. In higher organisms, the metabolism of individual tissues and organs is coordinated for the overall growth and well-being of the organism. Understanding the dependencies and rationale for multicellular metabolism is far from trivial. Here, we have advanced the use of AraGEM (a genome-scale reconstruction of Arabidopsis metabolism) in a multi-tissue context to understand how plants grow utilizing their leaf, stem and root systems across the day-night (diurnal) cycle. Six tissue compartments were created, each with their own distinct set of metabolic capabilities, and hence a reliance on other compartments for support. We used the multi-tissue framework to explore differences in the "division-of-labor" between the sources and sink tissues in response to: (a) the energy demand for the translocation of C and N species in between tissues; and (b) the use of two distinct nitrogen sources (NO(-) 3 or NH(+) 4). The "division-of-labor" between compartments was investigated using a minimum energy (photon) objective function. Random sampling of the solution space was used to explore the flux distributions under different scenarios as well as to identify highly coupled reaction sets in different tissues and organelles. Efficient identification of these sets was achieved by casting this problem as a maximum clique enumeration problem. The framework also enabled assessing the impact of energetic constraints in resource (redox and ATP) allocation between leaf, stem, and root tissues required for efficient carbon and nitrogen assimilation, including the diurnal cycle constraint forcing the plant to set aside resources during the day and defer metabolic processes that are more efficiently performed at night. This study is a first step toward autonomous modeling of whole plant metabolism.